{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Eigenvalue problem\n",
    "=================="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Matrix transformaion is written as\n",
    "\n",
    "$$ \\mathbf{y} = A\\mathbf{x} $$\n",
    "\n",
    "A different vector in the same direction can be written as scalar multiplication:\n",
    "\n",
    "$$\\mathbf{y} = \\lambda\\mathbf{x}$$ \n",
    "\n",
    "Equating these $y$s yields:\n",
    "\n",
    "$$ A\\mathbf{x} = \\lambda \\mathbf{x} \\Rightarrow (A - \\lambda I) \\mathbf{x} = 0$$\n",
    "\n",
    "$$\\det(A - \\lambda I) = 0 $$\n",
    "\n",
    "The eigenvalue problem can also be collected with $\\Lambda$ being a diagonal matrix containing all the eigenvalues and $X$ containing the eigenvectors stacked column-wise. This leads to the eigenvalue decomposition:\n",
    "\n",
    "$$  A X = X \\Lambda \\Rightarrow A = X \\Lambda X^{-1}$$\n",
    "\n",
    "with\n",
    "\n",
    "$$ \\Lambda = diag(\\lambda_i)$$\n",
    "\n",
    "If we try to find a similar decomposition with different constraints, we can write\n",
    "\n",
    "$$ A = U D V^{H} $$\n",
    "\n",
    "If $D$ is a diagonal matrix and $U$ and $V$ are [unitary](http://en.wikipedia.org/wiki/Unitary_matrix), this is the singular value decomposition.\n",
    "\n",
    "In Skogestad \n",
    "\n",
    "$$ A = U \\Sigma V^{H} $$\n",
    "\n",
    "$$ \\Sigma = diag(\\sigma_i)$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ipywidgets import interact"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotvector(x, color='blue'):\n",
    "    plt.plot([0, x[0,0]], [0, x[1,0]], color=color)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.patches as patches"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's investigate the properties of this matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = numpy.matrix([[4, 3],\n",
    "                  [2, 1]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The eigenvectors and eigenvalues can be calculated as follows. We also calculate the output vectors associated with a unit vector input in the eigenvector directions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[4, 3],\n",
       "        [2, 1]])"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "v = numpy.asmatrix(numpy.random.random(2)).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[0.89474813],\n",
       "        [0.44657115]])"
      ]
     },
     "execution_count": 80,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v = A*v\n",
    "\n",
    "v = v/numpy.linalg.norm(v)\n",
    "v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "lambdas, eigvectors = numpy.linalg.eig(A)\n",
    "ev1 = lambdas[0]*eigvectors[:, 0]\n",
    "ev2 = lambdas[1]*eigvectors[:, 1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The singular values determine the main axes of the translation ellipse of the matrix. Note that the `numpy.linalg.svd` function returns the conjugate transpose of the input direction matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "U, S, VH = numpy.linalg.svd(A)\n",
    "V = VH.H\n",
    "ellipseangle = numpy.rad2deg(numpy.angle(complex(*U[:, 0])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "412b4de2897149439ff0ce2a158eb4df",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>interactive</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "interactive(children=(FloatSlider(value=5.5, description='scale', max=10.0, min=1.0), FloatSlider(value=3.141592653589793, description='theta', max=6.283185307179586), Output()), _dom_classes=('widget-interact',))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.interactive>"
      ]
     },
     "execution_count": 83,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def interactive(scale, theta):\n",
    "    x = numpy.matrix([[numpy.cos(theta)], [numpy.sin(theta)]])\n",
    "    y = A*x\n",
    "\n",
    "    plotvector(x)\n",
    "    plotvector(y, color='red')\n",
    "    plotvector(ev1, 'green')\n",
    "    plotvector(ev2, 'green')\n",
    "    plotvector(V[:, 0], 'magenta')\n",
    "    plotvector(V[:, 1], 'magenta')\n",
    "    plt.gca().add_artist(patches.Circle([0, 0], 1, \n",
    "                                        color='blue', \n",
    "                                        alpha=0.1))\n",
    "    plt.gca().add_artist(patches.Ellipse([0, 0], S[0]*2, S[1]*2, \n",
    "                                         ellipseangle,\n",
    "                                         color='red',\n",
    "                                         alpha=0.1))\n",
    "    plt.axis([-scale, scale, -scale, scale])\n",
    "    plt.axes().set_aspect('equal')\n",
    "    plt.show()\n",
    "interact(interactive, scale=(1., 10), theta=(0., numpy.pi*2))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  },
  "widgets": {
   "state": {
    "89b320a62be44ebfa2ca3bd3e8610688": {
     "views": [
      {
       "cell_index": 13
      }
     ]
    }
   },
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}